# SI Figure 4
rm(list=ls())
gc()

# Temperature means and variances for which to estimate damages
mu <- seq(2,8,2)
sigma <- seq(0,0.5,0.01)

# Economic assumptions
eta = 1.45 # Marginal utility of consumption
gdp = 80000 # Aggregate consumption (billion USD)
pop = 7.5 # Population (billion people)
gdp.pc <- (gdp/pop)/1000 # Undistured per capita consumption (thousand USD).

# Generate random sample of temperature deviations
set.seed(1)
sample=100000000
sample=1000000
temp <- rnorm(sample,mean=0, sd=1)

# List all mean-variance combinations
d <- data.frame(mu=rep(mu,length(sigma)), sigma=rep(sigma,each=length(mu)))

# Compute expected damages for each mean-variance combination
library(doParallel)
cl <- makeCluster(2)
registerDoParallel(cl)
figure1 <- foreach(i=1:nrow(d), .combine='rbind') %dopar% {
	Temp <- temp*d$sigma[i] + d$mu[i]
	
	## Expected utility of consumption
	ExpD_Weitzman <- mean(((gdp.pc*(1/(1+(Temp/20.46)^2 + (Temp/6.081)^6)))^(1-eta) - 1)/(1-eta))
	c(d$mu[i],d$sigma[i],ExpD_Weitzman)
}
stopCluster(cl)
d <- as.data.frame(figure1,stringsAsFactors=F)
names(d) <- c("mu","sigma","ExpD_Weitzman")

UofC <- ((gdp.pc)^(1-eta) - 1)/(1-eta) # Utility of undisturbed consumption
scalar <- UofC-(((gdp.pc-0.001)^(1-eta) - 1)/(1-eta)) # Utility value of the marginal dollar

quartz(width = 5, height = 5, type = "pdf", "SI figure variance", file = "SI_figure_variance.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,1,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1), 1, byrow = TRUE), widths=c(3), heights=c(3))

plot(d$sigma,seq(0,4000000000000,length.out=nrow(d)),type="n",bty="L",xaxs="i",yaxs="i",ylab="",xlab="",xaxt='n',yaxt='n')
for (i in mu) {
	lines(d$sigma[d$mu==i],(pop*1000000000)*((UofC-d$ExpD_Weitzman[d$mu==i]) - (UofC-d$ExpD_Weitzman[d$mu==i & d$sigma==0]))/scalar,type="l", lty=1) # Economic loss in dollars
}
axis(side=1, at=c(0,.1,.2,.3,.4,.5), label=c("0",".1",".2",".3",".4",".5"), tick = TRUE,tck=c(.02))
axis(side=2, at=seq(0,10000000000000,1000000000000), label=seq(0,10,1), tick = TRUE,tck=c(.02), las=2)
mtext(side=3, "Trillion USD", outer=T, line=-4, cex=.9, at=0.1)

mtext(expression(paste("E(", Delta, "T)=2", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = (pop*1000000000)*((UofC-d$ExpD_Weitzman[d$mu==2 & d$sigma==.5]) - (UofC-d$ExpD_Weitzman[d$mu==2 & d$sigma==0]))/scalar,  adj = NA, padj = NA, cex = NA, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=4", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = (pop*1000000000)*((UofC-d$ExpD_Weitzman[d$mu==4 & d$sigma==.5]) - (UofC-d$ExpD_Weitzman[d$mu==4 & d$sigma==0]))/scalar,  adj = NA, padj = NA, cex = NA, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=6", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = (pop*1000000000)*((UofC-d$ExpD_Weitzman[d$mu==6 & d$sigma==.5]) - (UofC-d$ExpD_Weitzman[d$mu==6 & d$sigma==0]))/scalar,  adj = NA, padj = NA, cex = NA, col = NA, font = NA)
mtext(expression(paste("E(", Delta, "T)=8", sep="")), side = 4, las=2, line = 0, outer = FALSE, at = (pop*1000000000)*((UofC-d$ExpD_Weitzman[d$mu==8 & d$sigma==.5]) - (UofC-d$ExpD_Weitzman[d$mu==8 & d$sigma==0]))/scalar,  adj = NA, padj = NA, cex = NA, col = NA, font = NA)

mtext(expression(paste("Standard Deviation of ", Delta, "T", sep="")),side=1,outer=T, line=2.5)
dev.off()